Predicting Neutron Production from Cosmic-ray Muons 

Y-F. Wang, V. Balic, G. Gratta 
Physics Department, Stanford University, Stanford CA 94305, USA 

A. Fasso, S. Roesler 
Stanford Linear Accelerator Center, Stanford CA 94309, USA 

A. Ferrari 

INFN, Sezione di Milano, Via Celoria 16, 1-20133, Milano, Italy 
(February 4, 2008) 



O 
O 

(N 
C 

(N 



Fast neutrons from cosmic-ray muons are an important background to underground low energy 
experiments. The estimate of such background is often hampered by the difficulty of measuring 
and calculating neutron production with sufficient accuracy. Indeed substantial disagreement exists 
between the different analytical calculations performed so far, while data reported by different 
experiments is not always consistent. We discuss a new unified approach to estimate the neutron 
yield, the energy spectrum, the multiplicity and the angular distribution from cosmic muons using 
the Monte Carlo simulation package FLUKA and show that it gives a good description of most of 
the existing measurements once the appropriate corrections have been applied. 

PACS numbers: 96.40.Tv, 14.60. Pq, 25.20.-x, 25.40.Sc 
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I. INTRODUCTION 

Fast neutrons from cosmic-ray muons represent an im- 
portant background for low-energy underground experi- 
ments such as searches for proton decay and dark mat- 
ter, and low-energy neutrino oscillation. Unlike charged 
hadrons which can be tagged by a veto detector system, 
neutrons usually cannot be identified until they are fi- 
nally captured, mimicking the signal. The occasional 
neutron scattering with the sensitive material of the de- 
tector and the long lifetime of neutrons in the detector 
and surrounding materials further complicate the situ- 
ation. For example, the Palo Verde reactor neutrino 
oscillation experiment found such neutrons to be their 
dominant background Q], and a similar situation is ex- 
pected at the ultra-long baseline detector KamLAND || . 
Low-energy solar neutrino experiments such SNO and 
Borexino also have to estimate such backgrounds and the 
understanding of neutron backgrounds may be relevant 
in resolving the controversy between the CDMS || and 
DAM A jl| results on dark matter searches. Finally, low 
energy accelerator experiments at shallow depths, such as 
Karmen, LSND and OrLaND have also similar problems. 

Although the total neutron yield from cosmic muon 
spallation has been measured by several experiments, 
contradictory results are given in the literature 
Theoretical calculations [p|-pl| are also not consistent 
with each other and with data. The fact that primary 
neutrons, pions and protons can all produce secondary 
neutrons through hadronic interactions makes analytical 
calculations very difficult. A simple cascade model ]Tl| ] 
suggests that the number of nuclear cascade products 
such as neutrons, pions and protons, increases with the 



average muon energy approximately like E^' 7 . This for- 
mula agrees with most measurements ^-jj except the 
recent one by the LVD collaboration || . 

In addition to these problems with neutron yields, 
the few measurements of the neutron energy spec- 
trum [||,[l2 are not well reproduced by theoretical cal- 
culations [L|[l4|]. The interpretation of experimental 
data is complicated by the fact that the neutron en- 
ergy spectrum depends upon the muon spectrum that, 
in turn, is a function of the depth at which the mea- 
surement is carried-on and the geometrical configuration 
of the underground site. There is a broad range of re- 
sults reported in the literature. Barton [|l3| suggests 
that the spectrum of neutrons from hadronic cascade fol- 
lows E -1 / 2 between 10-50 MeV, while the spectrum of 
neutrons from n~ capture follows a flat spectrum up to 
100 MeV. Perkins Jl4| suggests that the neutron spec- 
trum from muon spallation follows E~ 16 . The com- 
bination of (9.7E- 1 / 2 + 6.0e~ B / 10 ) has been used in 
a measurement at a shallow site. It has also been 
suggested [[15| to use proton and neutron spectra follow- 
ing E~ 1S6 as measured at accelerators for photo-nuclear 
interactions p6| . Experimentally the Karmen experi- 
ment reported a visible energy spectrum following e ~ £ / 39 
for spallation neutrons Jl2j and the LVD experiment re- 
ported a visible energy spectrum following 

In order to put some order in this area, we have stud- 
ied the neutron yield, the neutron energy spectrum, the 
multiplicity and the angular distribution using the 1999 
version of the FLUKA JlTj Monte Carlo program which 
is expected to give an accurate description of all the pro- 
cesses involved. It is our intention to provide an unified 
approach and obtain a reliable estimate of neutron back- 
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ground for experiments at all depths. While we perform 
the simulation in a particular detector configuration, re- 
sults should be applicable to a range of substances, in- 
cluding rock and other shielding materials. 



II. MUON SPALLATION MODELS 

Fast neutrons from cosmic-ray muons are produced in 
the following processes: 

a) Muon interactions with nuclei via a virtual photon 
producing a nuclear disintegration. This process is 
usually referred to as "muon spallation" and is the 
main source of theoretical uncertainty. 

b) Muon elastic scattering with neutrons bound in nu- 
clei. 

c) Photo-nuclear reactions associated with electro- 
magnetic showers generated by muons. 

d) Secondary neutron production following any of the 
above processes. 

Processes b) and c) are reasonably well understood 
while a) and d) are the root of the difficulties described 
in previous calculations. Neutrons can be also produced 
from muons which stop and are captured, resulting in 
highly excited isotopes emitting one or more neutrons. 
This process is reasonably well understood and its con- 
tribution to total neutron yield can be calculated. All 
the experimental results referred to in this paper do not 
include these neutrons since they can be easily identified 
and eliminated. Neutrons produced from neutrinos are 
negligibly small at the depths considered, and thus are 
not discussed in this paper. 

The muon spallation process is schematically illus- 
trated in Figure 0. The desired [i — N cross section is 
then calculated as 



-N = 



du 



(1) 



where v = E — E', E and E' are energies of initial and 
final muons, and N 7 (z/) the virtual photon energy spec- 
trum. Theoretical calculations often treat the virtual 
photons according to the Weizsacker- Williams approxi- 
mation in which the passage of a charged particle 
in a slab of material produces the same effects as of a 
beam of quasi-real photons. A general expression of the 
Weizsacker- Williams formula is given in ref. fic|: 
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FIG. 1. The Feynman diagram of a muon spallation pro- 
cess. 



where m is the muon mass, and P and P' are momentum 
of initial and final muons. 

Since in the above approximation, it is assumed that 
the 7 — N cross section is the same for real and vir- 
tual photons, the measured 7 — N cross section can 
be used to calculate the fi — N cross section in Eq.(l). 
At low muon energy the situation is more complicated. 
Here, the virtuality of the photon becomes comparable 
to its energy and cannot be neglected. It follows that 
the Weizsacker- Williams approximation can no longer be 
used. In addition, the interaction of the virtual photon 
with the nucleus is a collective excitation of the nucleus 
(Giant Dipole Resonance, GDR) rather than a single 
photon-nucleon interaction. This implies that the GDR 
model would have to be applied to virtual photons intro- 
ducing further theoretical and technical complications. 
However, it might be reasonable to assume that neu- 
tron production by low-energy muon interactions is small 
as compared to neutron photoproduction by low-energy 
bremsstrahlung photons and adds therefore only a minor 
contribution to the total neutron yield. 

In addtional to these assumptions, there are a num- 
ber of problems associated with analytical calculations: 
first, they cannot reliably calculate all daughter products 
for every nucleus if the 7 — N interaction is very vio- 
lent so that the nucleus becomes highly excited; second, 
they cannot properly take into account secondary neu- 
tron production. Hence, while these calculations provide 
useful guidance and, at shallow depths, where hadronic 
shower effects are small, they may even give quantita- 
tively sound predictions p0| , in general they cannot be 
considered particularly reliable. 

Monte Carlo approaches are commonly used to prop- 
erly model hadronic cascades. Currently the most com- 
plete code to describe both hadronic and electromagnetic 
interactions up to 20 TeV is FLUKA 0. In this pro- 
gram, different physical models, or event generators, are 
responsible for the various aspects of particle production 
at different energies Jl9| . High-energy hadronic inter- 
actions arc described based on the Dual Parton Model 
followed by a pre-equilibrium-cascade model. In addi- 
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tion, models for nuclear evaporation, break-up of ex- 
cited fragments and 7-deexcitation treat the disintegra- 
tion of excited nuclei. Hadronic interactions of photons 
are simulated in detail from threshold (Giant Dipole Res- 
onance interactions) up to TeV-energies (Vector Meson 
Dominance Model). For nuclei up to copper, measured 
photo-nuclear cross sections in the low energy region are 
used |2l| . Hadronic interactions of muons are based on 
the Weizsacker- Williams approximation as formulated by 
Bezrukov and Bugaev p0| . A spectrum of virtual pho- 
tons is generated which interact with the nuclei similar to 
real photons. Due to the above theoretical and technical 
complications in the description of hadronic interactions 
of virtual photons at very low energies the simulation is 
restricted to photon energies above the delta resonance 
threshold. The implementation of hadronic interactions 
of muons has been shown to give reliable predictions for 
the MACRO experiment ||. 

In the following FLUKA is used to obtain a consis- 
tent and complete estimate of neutron production from 
cosmic muons. We model a simple cubic detector filled 
with liquid scintillator C n H2 n +2, where n is taken to be 
10. Muons with monochromatic energy are tracked in 
the detector and secondary neutrons and other hadrons 
are analyzed. 



III. NEUTRON YIELD 
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FIG. 2. Neutron production rate as a function of muon 
energy. The stars indicate the FLUKA simulation with a fit 
to the power law. The experimental points, with abscissa cor- 
responding to the average energy at the experiment's depth: 
A) 20 meter water equivalent (m.w.e.) j^), B) 25 m.w.e. p3[ , 
C) 32 m.w.e. by the Palo Verde experiment [0, D) 316 
m.w.e. ||, E) 750 m.w.e. |3, F) 3650 m.w.e. by the LVD 
experiment at Gran Sasso pf, and G) 5200 m.w.e. by the 
LSD detector at Mont Blanc H. 



The total neutron yield is probably the most measured 
quantity in our problem. There are many experimental 
results from different depths which can be compared with 
the model. 

All neutrons, either primary or secondary, are in- 
cluded. Double counting due to neutron scattering or 
neutron spallation is carefully avoided. Due to the lim- 
ited size of the detection volume in the simulation, some 
neutrons can escape, resulting in fewer secondary neu- 
trons. The total number of neutrons thus depends on 
the size of the detector. The detector size has to be lim- 
ited so that the muon energy-loss is small and the initial 
muon energy can be used as a constant. To correct this 
problem, we run the simulation with different detector 
sizes at each muon energy, and fit the neutron yield as a 
function of percentage of neutrons which escape the de- 
tector. An exponential behavior is found and the total 
number of neutrons can be extrapolated. The corrected 
neutron yield as a function of the muon energy is shown 
in Fig. ^. The neutron yield per muon can be fit as 

N„ = 4.14 x E°- 74 x 1CT 6 neutron/(/igcrrr 2 ) (3) 

where E is in GeV. This relation is consistent with E" 7 
law suggested in ref. JlT[ . 

While many experiments report their results as a func- 
tion of the detector's depth underground, clearly the 



proper physical parameter is the mean muon energy at 
the detector. The conversion between the two quanti- 
ties is not entirely trivial as the average energy depends 
upon the geometry of the overburden, particularly in the 
case of deep sites. Here for consistency we report all 
experimental results as a function of the average muon 
energy E M . For measurements performed at deep labora- 
tories such as Gran Sasso or Mont Blanc the conversion 
is given in the original papers ]^,|| , while for shallow sites 
we estimate the average muon energy-loss using the sim- 
ple relation |25f| and assuming a flat geometry: 

where E M is in GeV, a = (1.9 + 0.08 log(E M /m M )) x 
10~ 3 GeV/(gcm -2 ) and £ = 2.5 x 10 5 g/cm 2 in rock. 
All known experimental neutron-production rates are re- 
ported in Figure |^ along with the FLUKA predictions. 
The agreement is substantially better than obtained with 
previous calculations |pj^o|. 

Given this agreement we can proceed to extract from 
the Monte Carlo simulation information on the origin of 
the neutrons. This is of great interest as it can pro- 
vide hints for a better theoretical understanding of the 
processes involved. In Figure || we show the fractional 
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composition of the neutrons by origin. As expected the 
fraction of neutrons from primary processes such as muon 
spallation and real photo-nuclear interactions decreases 
with energy, while secondary neutrons from neutron and 
proton spallation, from pion absorption and from other 
minor processes such as A and £ decays, increase in rel- 
ative importance with energy. 



10 



O 



o 

TO 




FIG. 3. Origin of neutrons: a) direct muon spallation, b) 
real photo-nuclear disintegration, c) neutron spallation, d) 
proton spallation, e) n + spallation, f) tv~ spallation and cap- 
ture, g) others. 



The 7r + yield in muon spallation can be estimated from 
analytical calculations more reliably than that of neu- 
trons. This is because 1) secondary processes for neu- 
trons are more important; 2) there is an uncertainty in 
the direct neutron production from virtual or real pi- 
ons. Table | shows the comparison of our simulation 
with an analytical calculation JT(| which does not in- 
clude secondary pions from hadronic showers. In the 
experimental data (?]], 7r + 's are identified through their 
7T + — * /i + — > e + decays, and hadronic shower effects are 
not corrected for. In the pion case, good agreement be- 
tween data, FLUKA and the analytical calculation can 
be seen at shallow depth, where secondary hadrons from 
showers are not important. We obtained the pion yield 
per muon as: 



N 7r+ = 4.45 x E"- MU x 1Q-' pion/(^gcm^) 



(5) 



where E n 



is in GeV. The law appears to be univer- 



sal as suggested by Ryazhskaya [[□] and consistent with 
results from ref. 1241 . 



Depth(m) 
Energy(GeV) 


20 
10.3 


100 
22.4 


500 
80.0 


yield 


n 


7T + 


n 


7T + 


n 


7T+ 


anal. calc. juj 
FLUKA 
Data M 


0.87 
2.5 
3.0± 0.5 


0.30 
0.31 
0.35±0.07 


1.21 
3.9 


0.45 
0.52 


2.08 
11.0 


0.86 
1.51 



TABLE I. Neutron and n + yields (in units of 
10~ 5 /(gcm~ )) per muon at different depths. Note that 
the analytical calculation ]l(J does not include real pho- 
ton-nuclear disintegration and secondary particles. 



IV. NEUTRON ENERGY SPECTRUM 

The neutron energy spectrum is particularly controver- 
sial, with a wide range of results reported in theoretical 
calculations and in the few experimental measurements. 
In Figure ^ we show some of the energy spectra obtained 
with FLUKA. Each histogram is fitted to the universal 
empirical function 
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where A is a normalization factor, and 

B(E M ) = 0.52 - 0.58e-°- 0099£ " 



(6) 



(7) 



This simple function reproduces fairly well the FLUKA 
distributions with % 2 per degree of freedom of 3.9, 4.5, 
9.6 and 3.7 for 11, 20, 90 and 270 GeV respectively. 

We are aware of two direct measurements that can be 
compared to our calculations, as shown in Fig. [|. The 
Karmen collaboration measured neutron energy spec- 
trum up to 50 MeV Jl^l . Within this modest range, their 
parameterization N so fte _E,1 / 2 ' 1 +Nh a rde~ En / 39 , is in a rea- 
sonable agreement with our result. They attribute all the 
soft component to muon capture but it seems that muon 
spallation also produces some soft neutrons. The LVD 
experiment reported a E~ Y spectrum m up to 400 MeV, 
also in a reasonable agreement with our result. 

Some of the theoretical estimates of neutron energy 
spectrum are comparable to our results at low energies 
while our calculation gives results over a much wider 
range of energies. For example, the power laws E~ 1-6 |l4| ] 
and E~ lm O] agree with our results up to 400 MeV. 
Other functions, such as {9.7En 1/2 + 6.0e- £ ™/ 10 ) [§, or 
that suggested by Barton |0|, are significantly different. 



V. NEUTRON MULTIPLICITY AND ANGULAR 
DISTRIBUTION 

The neutron multiplicity is probably the least known 
quantity in the neutron production problem. In most 
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FIG. 4. Energy spectrum of neutrons at different muon energies together with our parameterization. We find \ 2 /NDF of 
3.9, 4.5, 9.6 and 3.7 for 11, 20, 90 and 270 GeV respectively. 
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FIG. 5. Comparison of measured neutron energy spectrum 
with our parameterization. 



yield. Similar results are also reported by Bezrukov et 
al. p3| ] . These numbers appear to be substantially larger 
than the FLUKA prediction of 2. While more data would 
be desirable to better understand this discrepancy it is 
possible that the data is affected by the incomplete ef- 
ficiency matrix for four or more neutrons, whose contri- 
bution, previously assumed small, seems significant from 
our simulation. 

There are no experimental results on the neutron an- 
gular distribution with respect to muons. Experimental 
results [^6j on 7 + 12 C — > p + X have been used |l5| ] 
as a first approximation of the neutron angular distribu- 
tion. The distribution is expected to be forward-peaked, 
smoothed somewhat by the contribution of secondary 
neutrons. Our simulation is well parameterized by the 
angular distribution 

dN - A 
dcosfl ~~ (l-cos0)°- 6 + B(E,J ( ' 

where B(E fl ) = 0.699E-" 136 . Figure shows the cosfl 
distribution at the usual set of energies together with the 
above parameterization. The x 2 P er degree of freedom 
are 0.51, 0.51, 0.72 and 0.72 for 11, 20, 90 and 270 GeV 
respectively. 

VI. SUMMARY 



cases muon spallation only happens once and produces 
only a few primary hadrons. But these hadrons can 
shower and generate secondary hadrons, including neu- 
trons. Using our simulation we have found that in some 
cases the number of secondary neutrons exceeds 50. The 
average number of neutrons is about 3 for a 11 GeV 
muon, and it increases to about 7 at a muon energy of 
385 GeV. Figure ^ shows the neutron multiplicity dis- 
tributions at different muon energies from FLUKA, to- 
gether with the universal empirical parameterization: 

= A (e- A ^ M + B(E^e- c ^ M ) (8) 
where M is the multiplicity, 

A(E^) = 0.085 + 0.54e" 075Et ' (9) 
27 2 

= I + 7 , 2e lo, 76 ^ (10) 

C(E„) = 0.67 + 1.4e-°- 12E * (11) 

The x 2 ' s P er degree of freedom are 0.6, 2.0,1.5 and 1.5 
for 11, 20, 90 and 270 GeV respectively. 

There are only few experimental results || on multi- 
plicity. The Palo Verde experiment at the shallow depth 
corresponding to a mean muon energy of about 16.5 GeV, 
observed a two-neutrons to one-neutron ratio between 5 
to 10, depending on the assumption on three neutron 



We obtained for the first time a complete description 
of neutron production from cosmic-ray muons using the 
FLUKA Monte Carlo program. Results have been com- 
pared with existing data and some analytical calculations 
reported in literature. With a few exceptions, our results 
agree well with data, and our predictions cover an en- 
ergy range much broader than what has been discussed 
before. 

Analytical paramctcrizations have been obtained for 
neutron yield, energy spectrum, multiplicity and angu- 
lar distribution. These formulae can be used as a start- 
ing point to estimate neutron backgrounds from cosmic- 
ray muons when the muon energy spectrum is known. 
For detailed calculation of muon-induced neutron back- 
grounds in low-energy underground experiments, a com- 
plete Monte Carlo simulation is needed to account for the 
correlation among different variables and can be setup 
using FLUKA. Intermediate processes such as neutron 
scattering before neutron spallation are also automati- 
cally taken into account in this method that has great 
promise to improve the quality of predictions. 
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FIG. 6. Neutron multiplicity at different muon energies together with our parameterization. We find \ 2 /NDF of 0.6, 2.0,1.5 
and 1.5 for 11, 20, 90 and 270 GeV respectively. 
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